setwd("/Users/ben/Documents/Papers/Current/EvoPVA/Quoll model R/Model inference")

###For both popsize and sex ratio###
## see below for just popsize ##

load("Kept_simsbothhalfpcED.RData")

stdz<-function(vec){
	(vec-mean(vec))/sd(vec)	
}
# Transform posterior samples
#re-generate test stats
obs<-c(766, 818, 470, 300, 64/(18), 87/(55), 46/(2), 30/(4), 3228, 4820, 2590, 2890, 134/(8), 167/(24), 84/1, 103/5)
temp<-gold[,8:23]
temp<-sweep(temp,2,obs)
temp<-temp^2
ss.pop.pb<-apply(temp[,c(1:4)], 1, sum) 
	ss.pop.as<-apply(temp[,c(9:12)], 1, sum) 
	std.ss.pop<-stdz(ss.pop.pb)+stdz(ss.pop.as) #sum standardised deviations
	ss.sr.pb<-apply(temp[,c(5:8)], 1, sum)
	ss.sr.as<-apply(temp[,c(13:16)], 1, sum)
	std.ss.sr<-stdz(ss.sr.pb)+stdz(ss.sr.as)
wghts<-sqrt(std.ss.pop^2+std.ss.sr^2) #weights as distance from circle centre
wghts<-1/wghts #inverse weights
rm(temp, obs)

par(mfrow=c(2,1))
plot(std.ss.pop, gold[,"beta"])
plot(std.ss.sr, gold[,"beta"])
plot(std.ss.pop, gold[,"alpha"])
plot(std.ss.sr, gold[,"alpha"])



beta.co<-lm(gold[,"beta"]~std.ss.pop+std.ss.sr, weights=wghts)$coefficients[2:3]
tf.beta<-gold[,"beta"]-beta.co[1]*std.ss.pop-beta.co[2]*std.ss.sr
hist(gold[,"beta"])
hist(tf.beta)

alpha.co<-lm(gold[,"alpha"]~std.ss.pop+std.ss.sr)$coefficients[2:3]
tf.alpha<-gold[,"alpha"]-alpha.co[1]*std.ss.pop-alpha.co[2]*std.ss.sr
hist(gold[,"alpha"])
hist(tf.alpha)

fsurv1.co<-lm(gold[,"fsurv1"]~std.ss.pop+std.ss.sr)$coefficients[2:3]
tf.fsurv1<-gold[,"fsurv1"]-fsurv1.co[1]*std.ss.pop-fsurv1.co[2]*std.ss.sr
hist(gold[,"fsurv1"])
hist(tf.fsurv1)

fsurv2.co<-lm(gold[,"fsurv2"]~std.ss.pop+std.ss.sr)$coefficients[2:3]
tf.fsurv2<-gold[,"fsurv2"]-fsurv2.co[1]*std.ss.pop-fsurv2.co[2]*std.ss.sr
hist(gold[,"fsurv2"])
hist(tf.fsurv2)

msurv.co<-lm(gold[,"msurv"]~std.ss.pop+std.ss.sr)$coefficients[2:3]
tf.msurv<-gold[,"msurv"]-msurv.co[1]*std.ss.pop-msurv.co[2]*std.ss.sr
par(mfrow=c(2,1))
hist(gold[,"msurv"])
hist(tf.msurv)

fec.co<-lm(gold[,"fec"]~std.ss.pop+std.ss.sr)$coefficients[2:3]
tf.fec<-gold[,"fec"]-fec.co[1]*std.ss.pop-fec.co[2]*std.ss.sr
par(mfrow=c(2,1))
hist(gold[,"fec"])
hist(tf.fec)

probd.co<-lm(gold[,"prob.d"]~std.ss.pop+std.ss.sr)$coefficients[2:3]
tf.probd<-gold[,"prob.d"]-fec.co[1]*std.ss.pop-fec.co[2]*std.ss.sr
par(mfrow=c(2,1))
hist(gold[,"prob.d"])
hist(tf.probd)



posteriors<-cbind(tf.beta, tf.alpha, tf.fsurv1, tf.fsurv2, tf.msurv, tf.fec, tf.probd)

post.summ<-function(vctr){
	temp<-summary(vctr)[c(4, 3, 1, 6)]
	pc2.5<-quantile(vctr, 0.025)
	pc97.5<-quantile(vctr, 0.975)
	SD<-sd(vctr)
	temp<-c(temp, pc2.5, pc97.5, SD)
	names(temp)[length(temp)]<-"SD"
	temp<-temp[c(1, 7, 2:3, 5, 6, 4)]
	temp
}

out<-t(apply(posteriors, 2, post.summ))

save(out, file="Posteriors halfpc ED Summary.RData")



###popsize only###

setwd("/Users/ben/Documents/Papers/Current/EvoPVA/Quoll model R/Model inference")

load("Kept_sims_poponlyhalfpcED.RData")

stdz<-function(vec){
	(vec-mean(vec))/sd(vec)	
}
# Transform posterior samples
#re-generate test stats
obs<-c(766, 818, 470, 300, 64/(18), 87/(55), 46/(2), 30/(4), 3228, 4820, 2590, 2890, 134/(8), 167/(24), 84/1, 103/5)
temp<-gold[,8:23]
temp<-sweep(temp,2,obs)
temp<-temp^2
ss.pop.pb<-apply(temp[,c(1:4)], 1, sum) 
	ss.pop.as<-apply(temp[,c(9:12)], 1, sum) 
	std.ss.pop<-stdz(ss.pop.pb)+stdz(ss.pop.as) #sum standardised deviations
	
wghts<-sqrt(std.ss.pop^2) #weights as distance from circle centre
wghts<-1/wghts #inverse weights
rm(temp, obs)

par(mfrow=c(2,1))
plot(std.ss.pop, gold[,"beta"])
plot(std.ss.pop, gold[,"alpha"])




beta.co<-lm(gold[,"beta"]~std.ss.pop, weights=wghts)$coefficients[2]
tf.beta<-gold[,"beta"]-beta.co[1]*std.ss.pop
hist(gold[,"beta"])
hist(tf.beta)

alpha.co<-lm(gold[,"alpha"]~std.ss.pop)$coefficients[2]
tf.alpha<-gold[,"alpha"]-alpha.co[1]*std.ss.pop
hist(gold[,"alpha"])
hist(tf.alpha)

fsurv1.co<-lm(gold[,"fsurv1"]~std.ss.pop)$coefficients[2]
tf.fsurv1<-gold[,"fsurv1"]-fsurv1.co[1]*std.ss.pop
hist(gold[,"fsurv1"])
hist(tf.fsurv1)

fsurv2.co<-lm(gold[,"fsurv2"]~std.ss.pop)$coefficients[2]
tf.fsurv2<-gold[,"fsurv2"]-fsurv2.co[1]*std.ss.pop
hist(gold[,"fsurv2"])
hist(tf.fsurv2)

msurv.co<-lm(gold[,"msurv"]~std.ss.pop)$coefficients[2]
tf.msurv<-gold[,"msurv"]-msurv.co[1]*std.ss.pop
par(mfrow=c(2,1))
hist(gold[,"msurv"])
hist(tf.msurv)

fec.co<-lm(gold[,"fec"]~std.ss.pop)$coefficients[2]
tf.fec<-gold[,"fec"]-fec.co[1]*std.ss.pop
par(mfrow=c(2,1))
hist(gold[,"fec"])
hist(tf.fec)

probd.co<-lm(gold[,"prob.d"]~std.ss.pop)$coefficients[2]
tf.probd<-gold[,"prob.d"]-fec.co[1]*std.ss.pop
par(mfrow=c(2,1))
hist(gold[,"prob.d"])
hist(tf.probd)



posteriors<-cbind(tf.beta, tf.alpha, tf.fsurv1, tf.fsurv2, tf.msurv, tf.fec, tf.probd)

post.summ<-function(vctr){
	temp<-summary(vctr)[c(4, 3, 1, 6)]
	pc2.5<-quantile(vctr, 0.025)
	pc97.5<-quantile(vctr, 0.975)
	SD<-sd(vctr)
	temp<-c(temp, pc2.5, pc97.5, SD)
	names(temp)[length(temp)]<-"SD"
	temp<-temp[c(1, 7, 2:3, 5, 6, 4)]
	temp
}

out<-t(apply(posteriors, 2, post.summ))

save(out, file="Posteriors halfpc ED Summary.RData")

